home *** CD-ROM | disk | FTP | other *** search
/ EnigmA Amiga Run 1996 February / EnigmA AMIGA RUN 04 (1996)(G.R. Edizioni)(IT)[!][issue 1996-02][Skylink CD III].iso / earcd / midi / gfft.lha / gfft-2.03 / source / gfft-2.03-source.lha / cfft.c < prev    next >
C/C++ Source or Header  |  1996-01-02  |  13KB  |  388 lines
  1. /***************************************************************************
  2.  *          Copyright (C) 1994  Charles P. Peterson                  *
  3.  *         4007 Enchanted Sun, San Antonio, Texas 78244-1254             *
  4.  *              Email: Charles_P_Peterson@fcircus.sat.tx.us                *
  5.  *                                                                         *
  6.  *          This is free software with NO WARRANTY.                  *
  7.  *          See gfft.c, or run program itself, for details.              *
  8.  *              Support is available for a fee.                      *
  9.  ***************************************************************************
  10.  *
  11.  * Program:     gfft--General FFT analysis
  12.  * File:        cfft.c
  13.  * Purpose:     fast Fourier analysis function with complex i/o array
  14.  * Author:      Charles Peterson (CPP)
  15.  * History:     4-June-1993 CPP; Created.
  16.  *              9-March-1994 CPP; added SAVE_TRIG_OPT sections.
  17.  *              6-Feb-95 CPP (1.31); Progress Requester
  18.  * Comment:
  19.  *   Inspired by the algorithm expressed in four1.c, page 411-412 in
  20.  *   Numerical Recipes in C  (William H. Press, Brian P. Flannery, Saul A.
  21.  *   Teukolsky, William T. Vetterling.  Cambridge University Press, 1988.)
  22.  *   In turn, inspiration for their work came from N. Brenner of Lincoln
  23.  *   Laboratories, Danielson, Lanczos, J. W. Cooley, J. W. Tukey, and many
  24.  *   others, including Fourier, for whom the underlying analysis is named.
  25.  *   This, however, is an original implementation without any copied code.
  26.  *   I hope it is much easier to read and understand than its predecessors.
  27.  *   If nothing else, it is very verbose, and doesn't require an intuitive
  28.  *   understanding of all trigonometric identities.  Things got a little
  29.  *   worse when I added an option to save the previously used trig values,
  30.  *   but most people would consider that (in principle) worthwhile.
  31.  *
  32.  *   Single precision floating point, except for trigonometric values
  33.  *   which are computed in double precision.
  34.  */
  35.  
  36. #define PROGRESS_INDICATOR
  37. #define SAVE_TRIG_OPT
  38. /*
  39.  * Comment the above lines out for easiest use outside GFFT.
  40.  *
  41.  * You will lose the optional saving of trig values, but (see note below)
  42.  * that really doesn't make a big difference in performance.
  43.  *
  44.  * Otherwise, when using this code outside GFFT you will have to duplicate
  45.  * the GFFT memory allocation function, or modify this file accordingly
  46.  * (e.g. with an error return if the memory demanded cannot be allocated).
  47.  *
  48.  * Note: When SAVE_TRIG_OPT is defined, as in GFFT, there is a BOOLEAN
  49.  * global variable SaveMemory which disables the saving of trig values
  50.  * when set.
  51.  */
  52.  
  53. #include <math.h>
  54.  
  55. #ifdef _AMIGA
  56. #ifdef _M68881
  57. #include <m68881.h>
  58. #endif
  59. #endif
  60.  
  61. #include "complex.h"
  62.  
  63. #ifdef SAVE_TRIG_OPT
  64. #include "gfft.h"       /* gmalloc safe allocation (longjmps on error) */
  65. #include "settings.h"   /* SaveMemory option */
  66. #endif
  67.  
  68. #ifdef PROGRESS_INDICATOR
  69. #include "wbench.h"
  70. #define PROGRESS_INCREMENT 1024
  71. long Inner_Loop_Count = 0;
  72. extern double Percent_Per_Loop;
  73. extern double Initial_Progress;
  74. static long next_progress_count = 0;
  75. #endif
  76.  
  77. void cfft (Complex_float datac[], long n, int isign)
  78. {
  79. /* ARGUMENTS:
  80.  *
  81.  * datac       input/output Complex_float data array
  82.  *             transformed according to isign (see below)
  83.  *
  84.  * n            number of complex data elements
  85.  *              *** n MUST BE A POWER OF 2 ***
  86.  *              *** THIS IS NOT CHECKED FOR ***
  87.  *
  88.  * isign    if 1, datac is replaced with its discrete Fourier transform
  89.  *              if -1, datac is replaced with n times its INVERSE
  90.  *                discrete Fourier transform
  91.  */
  92.  
  93. #ifdef SAVE_TRIG_OPT
  94.     static Complex_float *wkvector = NULL;
  95.     static last_n = 0;
  96. #endif
  97.  
  98. #ifdef PROGRESS_INDICATOR
  99.     next_progress_count = 0;
  100. #endif
  101.  
  102.  
  103. /* 
  104.  * This function is divided into 2 or 3 main blocks to minimize the scope
  105.  * of variables involved and hopefully to encourage better register use.
  106.  */
  107.  
  108. /* In the first block we sort the input array of complex numbers into
  109.  * bit reversal order (e.g. the value at binary index ...001 is swapped
  110.  * with the value at binary index 100...)  
  111.  *
  112.  * While i counts forwards from 1, j counts "inversely" from binary 100...
  113.  *
  114.  * Skip i==0 and i==n-1 (They never get swapped anyway)
  115.  */
  116.  
  117.     {
  118.     long const halfn = n / 2;
  119.     long const nminus1 = n - 1;
  120.     long register j = halfn;
  121.     long register i = 1;
  122.  
  123.     for (; i < nminus1; i++)
  124.     {
  125.         if (j > i)    /* If move required, and not done before */
  126.         {
  127.         CF_SWAP (datac[i], datac[j]);
  128.         }
  129.     /* 
  130.          * The following sub-block is for inverse counting.
  131.          * To count inversely, we clear highest set bits until reaching
  132.          * the first clear bit, then set it.
  133.          */
  134.         {
  135.         long register scanbit = halfn;
  136.         while (j >= scanbit && scanbit > 0)
  137.         {
  138.             j -= scanbit;
  139.             scanbit >>= 1;
  140.         }
  141.             j += scanbit;
  142.  
  143.         }   /* End inverse counting sub-block */
  144.  
  145.     }   /* End for (i = 1, 2, ..., n-1) */
  146.  
  147.     }   /* End sorting block */
  148. /*
  149.  * The following block(s) are the Danielson-Lanczos section of the routine.
  150.  * Here we combine the transforms, initially of one value each, by two
  151.  * at a time giving us transforms twice as large, until everything has
  152.  * been combined into one big transform.
  153.  *
  154.  * Thanks to clever sorting, each transform of 'even' elements is
  155.  * followed by a transform of 'odd' elements during each combination.
  156.  *
  157.  * mathematical summary (in which = denotes equality or definition):
  158.  * (A transform of length one simply copies input to output)
  159.  * F[k] = F[k][even] + (w ** k) * F[k][odd]
  160.  * (w ** k) = - (w ** (k + N/2))
  161.  * w = e ** (2 * PI * i / N)
  162.  * e ** (i * theta) = cos (theta) + i * sin (theta)
  163.  * i = sqr (-1)
  164.  * e and PI are the named transcendental numbers
  165.  * cos and sin are the named trigonometric functions
  166.  * N is the number of complex samples; k is a number 0..N-1
  167.  */
  168.  
  169. #ifdef SAVE_TRIG_OPT
  170.     if (SaveMemory)
  171. #endif
  172.  
  173.    /*
  174.     * This is the unaccelerated version
  175.     * The trig constants are computed as they are needed
  176.     * each time cfft is called.
  177.     */
  178.  
  179.     {
  180.     long osize;    /* Size of old transforms being combined */
  181.     long nsize;    /* Size of new transforms being produced */
  182.     Complex_double w;
  183.     Complex_double wtemp;
  184.     Complex_double wk;        /* (w ** k) */
  185.     Complex_float wk_times_f_k_odd;
  186.     Complex_float wkfloat;    /* Floating copy for innermost loop */
  187.  
  188.     double theta;
  189.     long ieven, iodd, k;
  190.  
  191.     for (osize = 1,nsize = 2;  osize < n;  osize = nsize,nsize *= 2)
  192.     {
  193.         theta = 2.0 * PI / (isign * nsize);  /* PI is from math.h */
  194.         w.real = cos (theta);
  195.         w.imaginary = sin (theta);
  196.         wk.real = 1.0;
  197.         wk.imaginary = 0.0;
  198.         wkfloat.real = 1.0;
  199.         wkfloat.imaginary = 0.0;
  200.  
  201.         for (k = 0; k < osize; k++)
  202.         {
  203.         for (ieven = k; ieven < n; ieven += nsize)
  204.         {
  205.             iodd = ieven + osize;
  206.         /*
  207.                  * Note that the the upper and lower parts of each new
  208.                  * transform being produced use the same components--
  209.                  * taken through a 'second cycle' since they are periodic
  210.          * in their original N (nsize/2).  Only W**k differs,
  211.          * and, fortuitously, W**k = -W**(k+nsize/2).  Also
  212.          * fortuitously, the input slots for F_k_even and F_k_odd
  213.          * are separated by nsize/2, as are the output slots for 
  214.                  * F_k and F_(k+nsize/2), so we can use and write over the 
  215.                  * same slots during each pass in the inner loop, using the
  216.                  * same value of W**k (we subtract instead of negating).
  217.          */
  218.             C_MULT (wkfloat, datac[iodd], wk_times_f_k_odd);
  219.             C_SUB (datac[ieven], wk_times_f_k_odd, datac[iodd]);
  220.             C_ADD (datac[ieven], wk_times_f_k_odd, datac[ieven]);
  221.  
  222.         } /* end for (sums using same wk factor) */
  223.  
  224.         /*  C_MULT (wk, w, wk) but must use temp or i/o coincide */
  225.  
  226.         C_MULT (wk, w, wtemp);
  227.         wkfloat.real = (float) (wk.real = wtemp.real);
  228.         wkfloat.imaginary = (float) (wk.imaginary = 
  229.                          wtemp.imaginary);
  230.  
  231.         } /* end for (k = 1, 2, ..., osize) */
  232.  
  233.     } /* end for (osize, nsize...) */
  234.  
  235.     } /* end of Danielson-Lanczos section (SaveMemory slow version) */
  236.  
  237. #ifdef SAVE_TRIG_OPT
  238.     else /* if (!SaveMemory) */
  239.     {
  240.  
  241. /* To speed things up under most cases (where the same N is used for
  242.  * more than one transform, we are going to first create the vector of
  243.  * trigonometrically derived wk constants.  If N is unchanged, and the
  244.  * the vector has already been created, we simply re-use the previous
  245.  * vector.  If it weren't for this section (which uses GFFT allocation
  246.  * routines) there would be no memory allocated in cfft, and no need for
  247.  * the gfft.h header file).  Note that this also increases the dynamic
  248.  * memory requirements of gfft significantly because the vector necessary
  249.  * must be large enough for 2*n complex floating point numbers, making it
  250.  * one of the biggest single dynamic memory uses--particularly for large N
  251.  * which some people might be interested in.  So, I have decided
  252.  * to make this feature optional, if included at all.  As currently
  253.  * implemented, it also requires the 'safe' gmalloc, which longjmps
  254.  * on error (so as not to confuse the interface here and elsewhere).
  255.  *
  256.  * 11-March-94 CPP; I'm somewhat disappointed with the 5-7% overall
  257.  * performance improvement attributable to this modification, in spite of
  258.  * the fact that lstat shows cfft as using 67% (58% accelerated) of the
  259.  * total time for a 30 sec batch spectrum analysis (run on a vanilla 68000
  260.  * with no FFP even).  But, I should have anticipated this minor change;
  261.  * the trigs become a relatively small part of the overall calculations.
  262.  * Don't let anyone tell you any different.  (Meanwhile, the dreaded
  263.  * ok_spectrum and ok_read, which look so horribly inefficient, were only
  264.  * taking around 2% of the total time each.)  The number one line is the
  265.  * complex subtraction following the main fft complex multiplication; it
  266.  * probably waits there for the complex multiplication to finish.  That,
  267.  * not surprisingly, it the chief bottleneck.  As it should be.
  268.  *
  269.  */
  270.  
  271.     if (n > 1 && (!wkvector || n != last_n))
  272.     {
  273.         long osize;    /* Size of old transforms being combined */
  274.         long nsize;    /* Size of new transforms being produced */
  275.         Complex_double w;
  276.         Complex_double wk;          /* (w ** k) */
  277.         Complex_double wtemp;
  278.  
  279.         double theta;
  280.         long k;
  281.         long vsize;  /* vector size */
  282.         long ntemp;
  283.         long wki = 0;
  284.     /*
  285.      * The vector size is N + log2(N) - 1 where log2 is 'log base 2'
  286.      * we calculate this using integer arithmetic
  287.      */
  288.         vsize = -1;
  289.         for (ntemp = n; ntemp > 1; ntemp /= 2) vsize++; /* log base 2 */
  290.         vsize += n;
  291.  
  292.         if (wkvector)
  293.         {
  294.         gfree (wkvector);
  295.         }
  296.         wkvector = gmalloc (vsize * sizeof(Complex_float), 
  297.                 NOTHING_SPECIAL);
  298.  
  299.         for (osize = 1,nsize = 2;  osize < n;  osize = nsize,nsize *= 2)
  300.         {
  301.         theta = 2.0 * PI / (isign * nsize);  /* PI is from math.h */
  302.         w.real = cos (theta);
  303.         w.imaginary = sin (theta);
  304.  
  305.         wk.real = 1.0;
  306.         wk.imaginary = 0.0;
  307.  
  308.         wkvector[wki].real = 1.0;
  309.         wkvector[wki++].imaginary = 0.0;
  310.  
  311.         for (k = 0; k < osize; k++)
  312.         {
  313.         /*  C_MULT (wk, w, wk) but must use temp or i/o coincide */
  314.  
  315.             C_MULT (wk, w, wtemp);
  316.             wkvector[wki].real = (float) (wk.real = wtemp.real);
  317.             wkvector[wki++].imaginary = (float) (wk.imaginary = 
  318.                              wtemp.imaginary);
  319.  
  320.         } /* end for (k = 1, 2, ..., osize) */
  321.  
  322.         } /* end for (osize, nsize...) */
  323.  
  324.         last_n = n;
  325.  
  326.     } /* end of trig constants computation block */
  327.  
  328.     /*
  329.      * Now, we compute this fft
  330.      */
  331.         {
  332.         long osize;    /* Size of old transforms being combined */
  333.         long nsize;    /* Size of new transforms being produced */
  334.         Complex_float wk_times_f_k_odd;
  335.  
  336.         long ieven, iodd, k;
  337.         register Complex_float *wkp = wkvector;
  338.  
  339.         for (osize = 1,nsize = 2; osize < n;  osize = nsize,nsize *= 2)
  340.         {
  341.         for (k = 0; k < osize; k++)
  342.         {
  343.             for (ieven = k; ieven < n; ieven += nsize)
  344.             {
  345.             iodd = ieven + osize;
  346.         /*
  347.          * Note that the the upper and lower parts of each new
  348.          * transform being produced use the same components--
  349.          * taken through a 'second cycle' since they are periodic
  350.          * in their original N (nsize/2).  Only W**k differs,
  351.          * and, fortuitously, W**k = -W**(k+nsize/2).  Also
  352.          * fortuitously, the input slots for F_k_even and F_k_odd
  353.          * are separated by nsize/2, as are the output slots for 
  354.          * F_k and F_(k+nsize/2), so we can use and write over the 
  355.          * same slots during each pass in the inner loop, using the
  356.          * same value of W**k (we subtract instead of negating).
  357.          */
  358.             C_MULT (*wkp, datac[iodd], wk_times_f_k_odd);
  359.             C_SUB (datac[ieven],wk_times_f_k_odd, datac[iodd]);
  360.             C_ADD (datac[ieven],wk_times_f_k_odd, datac[ieven]);
  361.  
  362. #ifdef PROGRESS_INDICATOR
  363.             if (++Inner_Loop_Count >= next_progress_count)
  364.             {
  365.                 next_progress_count += PROGRESS_INCREMENT;
  366.                 progress_requester__update 
  367.                   ((int) (Initial_Progress + (Inner_Loop_Count *
  368.                    Percent_Per_Loop)));
  369.             }
  370. #endif
  371.             } /* end for (sums using same wk factor) */
  372.  
  373.             wkp++;
  374.  
  375.         } /* end for (k = 1, 2, ..., osize) */
  376.  
  377.         wkp++;
  378.  
  379.         } /* end for (osize, nsize...) */
  380.  
  381.     } /* end of transform computation */
  382.  
  383.     } /* end of !SaveMemory (fast memory hog) version */
  384.  
  385. #endif /* whether SAVE_TRIG_OPT defined or not */
  386.  
  387. } /* end of cfft */
  388.